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Abstract. In this letter, we demonstrate that oscillations and excitable behavior can 
be imparted to a chemical network by coupling the network to an active hydrogel. We 
discuss two mechanisms by which the mechanical response of the gel to the embedded 
chemical reactant provides feedback into the chemistry. These feedback mechanisms 
can be applied to control existing chemical oscillations as well as create new oscillations 
under some conditions. We analyze two model systems to demonstrate these two 
effects, respectively: a theoretical system that exhibits no excitability in the absence 
of a gel, and the Oregonator model of the Belousov-Zhabotinsky reaction in which 
the metal catalyst is intercalated into the polymer network. This work can aid in 
designing new materials that harness these feedbacks to create, control, and stabilize 
oscillatory and excitable chemical behavior in both oscillatory and non-oscillatory 
chemical networks. 


PACS numbers: 82.35.Jk, 82.70.Gg, 82.33.Ln, 05.45.a 

1. Introduction 

Nonlinear chemical networks are popular candidates for model systems of complex 
biological and biochemical phenomena. Such networks can harness a balance of positive 
and inhibitory feedback controls to generate oscillatory and excitable behavior in both 
time and space. The complex networks and idealized model systems serve for example 
as analogs for studying various biological processes such as clocks, embryo development, 
and nervous system behavior P1211111 El IE]. 

Recently, there is a growing interest in utilizing these networks to go beyond 
studying natural systems to generating new materials with smart and biomimetic 
properties piHiiniiin]. 

For example, it has been shown mm that by anchoring the metal catalyst of the 
Belousov-Zhabotinsky reaction to the polymer gel network of a hydrogel, the oscillatory 
redox reaction induces swelling and deswelling of the gel. This swelling can be utilized 
to induce biomimetic behaviors such as peristaltic motion 
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Figure 1. Schematic of chemomechanical and mechanochemical influence in a 
chemically-coupled hydrogel. The box indicates the bare chemical network. Solid and 
broken arrows indicate “dynamic” and “passive” influence, respectively, as described 
in the main text. 

These advances show that coupling a hydrogel to a chemical network induces new 
and useful behaviors in the gel. Herein we argue that the complementary effect is of equal 
interest na El na ca [H]: the mechanical gel contributes additional feedback to the 
chemical network. The additional mechanical feedback controls the stability of chemical 
limit cycles, and in some cases, introduces new excitable and oscillatory behavior that 
is otherwise impossible without mechanical coupling. Such behavior has been predicted 
and demonstrated mm using a bistable reaction driving a pH sensitive gel. In those 
experiments, the swelling and deswelling of the gel modulates the mixing rate between 
the gel pore-space and the well-mixed surrounding reservoir. The chemistry is thereby 
coupled to the gel via transport and mixing processes. 

The above work exploring transport induced chemical oscillations, as well as work 
on chemically coupled BZ-gels, serve as inspiration for the work in this letter. Here, 
we demonstrate the roles of mechanical feedback in two model chemical networks in 
which the chemistry is directly coupled to the polymer network itself. First, we analyze 
a simple theoretical network that exhibits excitable and oscillatory behavior only when 
mechanically and chemically coupled to the surrounding medium. Second, we show that 
stability and magnitude of preexisting characteristic oscillations can be controlled by 
varying the material parameters of the coupled hydrogel. 

We then discuss the several mechanisms by which the swelling and deswelling of 
the polymer network contribute feedback. By coupling a chemical network to an active 
hydrogel, the additional feedback mechanisms can be taken advantage of in designing 
systems with tightly controlled constant, excitable, and oscillatory behaviors. 

2. Mechanically coupled chemical network 

To understand the dynamical coupling described above, we consider a chemical reaction 
ongoing within the pore space of a hydrogel. The reaction is characterized by two 
concentrations, u and v, both measured with respect to the total volume (solvent and 
polymer), rather than solvent volume [IHj. The condition for this approximation is 
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that the microscopic mixing lengthscale of polymer molecules is much smaller than the 
distance between reacting molecules. The chemical network is indicated in Figure [U by 
the box and is characterized by two kinetic equations: 

^ = F(n,v) 

The chemical reaction is coupled to the enveloping hydrogel by embedding one of 
the reactants into the polymer network itself. For example, in BZ-coupled-hydrogels [TT] . 
the metal catalyst is anchored via a bipyridine linker. As the concentration of coupled 
reactant changes the polymer network becomes more or less hydrophobic [H]. The 
internal stress, cr, therefore depends on the concentration of polymer-bound reactant, v. 
In this way, the physical gel is influenced by the chemistry, as indicated by the arrow 
pointing from the box to a in Figure [H 

Specifically, a depends on both v and the polymer volume fraction, 0. We model 
the dependence with the Flory-Huggins osmotic pressure for an isotropic elastic gel, 
modified to account for the interaction between solvent and anchored reactant|18|: 

^ =-(0 + log(l-0)+ x(T)(/) 2)+ xW(/) + cono j, (2) 

where y is the temperature dependent solvent-polymer coupling parameter, cq is the 
density of cross-linkers, Vq is the volume of each statistical polymer subunit, 00 is a 
reference polymer density, and y* is the parameter measuring the coupling between 
solvent and anchored reactant. The gel responds to changes in a by swelling and 
deswelling. Under the approximation that the gel is homogeneous in that reactants 
and polymer are evenly distributed throughout the gel, we can write 

^ = (3) 

where M measures the mechanical susceptibility, or mobility, of the gel. For further 
discussion of a and derivation of Eq. (3), we refer the reader to the Supplementary 
Information. 

As the internal pressure varies, the gel expels and absorbs fluid thereby shrinking 
and swelling. This process contributes feedback from mechanics to the chemical 
kinetics in two ways. First, non-limiting reactants with concentrations fixed relative 
to the outside solvent are diluted by the varying polymer fraction, 0. Suppose the 
concentration of such a reactant is A in the bath. When measured within the gel relative 
to total volume, the concentration becomes (1 — 0)A [IH]. Therefore, each factor of A 
in the kinetics is replaced by the 0-dependent rate (1 — 0)A. 

Second, the concentration of polymer increases as the gel contracts, causing an 
increase in the concentration of bound reactant, v. Conversely, as the gel swells, the 
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Figure 2. Damped oscillatory behavior demonstrated by the theoretical network. The 
top panel gives the concentration of v both with the coupling (black line) and without 
(grey line). The former exhibits an initial jump whereas the latter rapidly decays to 
'T’ = 0, as shown in the detailed inset. The bottom panel gives the polymer density (j). 
The model parameters are a = 0.1, /3 = 0.2, M = 1, cqVq = 3 x 10“^, T = 306, and 
y* = 0.045. 


concentration v decreases. On the other hand, the concentration of dissolved reactant, 
decreases as the gel contracts because the polymer behaves as a dilutant. With this 
additional “dilution” effect, the time dependencies in Equation (1) become 


du 

dt 

dv 

dt 


F{u{t),v{t),(t>{t)) - 
G{u{t),v{t),(p{t)) + v^. 


(4) 


In the homogeneous approximation diffusion is neglected as well as dilution of the 
dissolved species at the edges of the gel due to mixing of solvent with the reservoir. 
These effects are known to induce oscillations in some cases, as discussed in [IHl [IS] and 

m. 

In addition, we show that oscillations can be created and controlled in the absence 
of mixing processes. We apply the above coupled gel model to two chemical networks. 
The first is a theoretical network of our design that exhibits no oscillatory or excitable 
behavior in the absence of an active gel: 


F{u, n, 0) = 1 — u{l — (j)Y 
G{u, v) = auv — (dv, 


where a and (3 are dimensionless constants. However, we have introduced chemo- 
mechanical feedback by assuming that the degradation rate of the dissolved component, 
u, depends on non-limiting reactants and therefore inherits factors of (1 — 0). As we 
will see, this additional coupling introduces excitable behavior to the system. 

The second model chemical network is a simplified form of the Oregonator, an 
approximate description of the Belousov-Zhabotinsky (BZ) reaction P . For a thorough 
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Table 1. Parametric dependence of coupled Oregonator 
Parameter Behavior under increase 

X* Switches from small-type to large-type oscillations and generally increases amplitude 

cqUo Decreases amplitudes. Switches from large-type to small-type oscillations 

T Removes oscillations (Low y*: decreases amplitudes; High y*: increases amplitudes) 

M Removes oscillations at low M, but introduces oscillations at high M 


description of the BZ reaction and Oregonator, we direct the reader to pun], and I2D1. 
In experimental studies, the chemical reaction is coupled to the hydrogel by linking 
the Ruthenium catalyst to the polymer m- Accounting for fixed concentrations of 
nonlimiting reagents, the kinetics for the reactants are given by |18| : 

F(u, «» = (!- 

eG{u, V, 0) = (1 — (f)Yu — (1 — (f))v, 

where /, q, and e are well known BZ specific parameters. 

The above dynamical systems are represented in Figure [T] which shows the direction 
of influence in the coupling among the chemistry (the box), the internal pressure a, and 
the physical gel state 0. There are two classes of coupling between the mechanics of 
the gel and the chemistry. The first, shown in solid arrows, we refer to as “dynamic” 
because dependencies are in the time derivatives of the variables. The second, shown 
in broken arrows, we refer to as “passive” because the dependencies are functional, but 
not explicitly in the time derivatives. The former are scaled by the mobility of the 
gel, M, whereas the latter are controlled by the polymer-catalyst coupling strength, 
X*. Thus, the model described above presents a chemo-mechanically and mechano- 
chemically coupled system. 

3. Results 

3.1. Theoretical network 

To show that coupling a non-oscillatory chemistry to a gel introduces new behaviors, 
we numerically integrate the theoretical model, Eqs. (4)-(5), for various values of a 
and 13. We find excitable behavior for many parameter sets, including that shown in 
Figure [2], where the chemistry responds to a perturbation from the stationary state 
with an initial jump followed by damped oscillations. In the absence of coupling (where 
(1 — 0) is constant), the chemical system quickly and smoothly returns to the steady 
state following a perturbation. Results of a linear stability analysis confirm these 
observations, that the uncoupled system will never oscillate and the coupled system can 
exhibit excitable jumps followed by damped oscillations. However, in the Supplementary 
Information, we present a similar three component model that exhibits undamped 
oscillatory behavior. Oscillatory behavior was previously reported for a non-oscillatory, 






Coupling chemical networks to hydrogels controls oscillatory behavior 


6 



Figure 3. Trajectories of system in v and cj) originating at two distinct initial 
conditions (green and black dots) that approach stable oscillations near stationary 
states (red dots). The green trajectory forms a simple limit cycle around the deswollen 
fixed point whereas the black trajectory forms a ringing oscillation around the swollen 
fixed point. The parameters are / = 0.6, e = 0.15, M = 12.6, cqVq = 3 x 10“^, 
T = 303, and y* = 0.0134. 

but bistable, chemical network when immersed within a hydrogel ini In that system, 
the oscillatory behavior relies upon diffusive and advective communication between the 
bath and a bistable chemistry. By coupling the reaction directly to the hydrogel, the 
requirements for both bistability and spatial resolution are therefore relaxed. 

3.2. Parametric control of oscillatory stability in Oregonator 

In addition to creating new chemical oscillations, mechanochemical feedbacks modify 
the stability of existing oscillatory behavior. 

Stationary states in the two-variable Oregonator are known to exhibit Hopf 
bifurcations under variations in the BZ parameters, leading to stable oscillations of 
u and V. In the case of the BZ-coupled gel, linear stability analysis demonstrates that 
under the assumption of large mobility, M, variations in y* can destabilize the stationary 
state [U]. 

We expand on that work by searching for stable oscillations while varying 4 
parameters describing our coupled system which can be controlled experimentally via 
the design, condition, and preparation of the gel and reaction. The parameters are 
the mobility, M, solvent-catalyst coupling coefficient, y*, crosslink density cqVq, and 
temperature T. For each set of parameters, we integrate Equation (4) from many initial 
conditions to find the magnitude of stable limit cycles. The general behaviors observed 
by varying the four parameters are summarized in Table [H 

As each parameter is increased, the stability and amplitude of oscillations change. 
These results are useful for designing BZ-coupled gels that exhibit engineered behaviors 
under varying conditions. More generally, the results demonstrate that the coupling 
causes the chemistry to adopt non-trivial and powerful dependencies on mechanical 
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Figure 4. Two types of mechano-chemical feedback (in red) modifying the three- 
variable Oregonator (in black). The first type depends on (p and the second depends 
on (1 — ()>)”. Reactants and influence lead into the reactions represented by black dots. 
The exponents appropriate for the BZ reaction are a = 2, b = d = 1, and c = 3. 

parameters. Therefore, the feedback mechanisms that are associated with these 
parameters offer new ways to control oscillatory chemical networks. 

The addition of the mechanical variable, 0, to the two-variable Oregonator permits 
new types of oscillatory behavior. If the mechanical dynamics vary over a timescale 
different from the chemical kinetics, ringing oscillations are observed as shown by the 
black curve in Figure |3l Each pulse consists of a slow compression of the gel followed by 
an excitation of the chemistry. The chemistry then rapidly oscillates around the u — v 
stationary point as (f) slowly decreases towards the global stationary point. 

The ringing behavior depends on initial conditions; given one set of parameters, 
multiple stable oscillations are available and can be selected by initializing the system 
in the swollen or deswollen state, as shown by black and green curves in Figure El 

4. Discussion 

Above, we describe how two novel feedback mechanisms are introduced by coupling 
gel and chemistry. The two mechanisms are illustrated in red in the context of the 
BZ reaction in Figure IH The species Br“ is the inhibitor in the Oregonator, which 
is assumed to be at rapid equilibrium in our two-variable model. The first feedback 
mechanism is the modification of the kinetics due to the effective dilution of nonlimiting 
reactants by polymer. This gives rise to the (1 — (pY terms. The second mechanism 
is the dynamic dilution and concentration of reactant concentrations due to the time 
derivative of cp. This effect is notated in Figure E] by arrows controlled by (p. 

The factors of (1 — (pY influence the kinetics both dynamically and at steady state. 
From one perspective, the quantity (1 — 0) acts as a “reactant” in its own right as 
reflected in Eq.(3). When v increases, (1 — 0) increases to approach osmotic equilibrium. 
It is included in the system of differential equations almost on the same footing as u and 
V. In this interpretation, the new “reactant” acts as a bridge, as shown schematically 
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Figure 5. Chemical schematic illustrating the role oi{l — (p) in the theoretical network. 


in Figure El 

From another perspective, the dilution of the non-limiting reactants, as captured 
by factors of (1 — (p)^, simply modify the kinetic rates. Therefore, the linear stability of 
the stationary states depends on the particular exponents. As the value v is perturbed 
from the stationary state, the equilibrium value for (1 — (p) is perturbed and the gel 
expands or contracts at a rate proportional to the mobility. The new value of (1 — 0) 
feeds back to various reactions with a strength related to the exponent, whose value 
is determined by the stoichiometry. This feedback then either damps or magnifies the 
initial perturbation. 

We can quantify this effect by computing the largest real part of the eigenvalues of 
the linearized system at stationary state as a function of the exponents and the mobility, 
M (see Figure El (top row)). To isolate the effect of feedback, we normalize factors of 
(1 — <p) in Equation (4) by the stationary value, (1 — 0ss). We find a nonmonotonic 
dependence on mobility where for small M the stationary state is stabilized for a > b, 
but the opposite is true for large M. This rich variety of feedback mechanisms introduces 
nonintuitive effects as the mechanical mobility is adjusted. 

We next examine Figure 0] to anticipate how the cp terms influence the overall 
dynamics. Consider for example the rising side of a pulse, in which the gel is oxidizing. 
As V increases, the gel swells and 0 < 0. The rate v consequently decreases. This 
contributes an “inertia” to the system, in that it resists change in the oxidation 
state. Stationary states are unaffected, but oscillations, excitations, and transients 
will experience a “drag”. This mechanism is responsible for the damped oscillations 
exhibited in Figure El The alternate model presented in the Supplementary Information 
swaps the sign of the drag term, thereby stabilizing the oscillations and avoiding the 
overdamped situation. Figure El (bottom row) shows how the stability of the stationary 
state in the Oregonator is controlled by these inertial terms in a nonmonotonic way. In 
designing a chemical network, this drag effect can be utilized to resist rapid changes in 
the reactant concentrations, thus contributing an inductance-like effect in the chemical 
network. 
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Figure 6. Stability of the system at stationary state for various mobilities and 
exponents a and b. Each square gives a and b = 0,1,2,3 (for all cases c = d = 0) 
for one value of the mobility. Blue indicates at least one eigenvalue has positive real 
part (unstable), whereas red indicates all negative real parts (stable). All factors of 
(1 — (/>)” are normalized. The broken white line separates positive from negative. Top 
and bottom rows correspond to systems lacking and containing the (j) term, respectively. 


5. Conclusions 

In this work we have furthered the understanding of how chemomechanical and 
mechanochemical couplings control and induce chemical oscillations in coupled 
hydrogels. It has been suggested that such systems offer applications as biomimetic 
devices and micro-actuators [H]. The additional control described here could be critical 
in engineering finely tuned chemical oscillations in such devices as well as ensuring 
behavioral robustness. Our analysis helps identify perturbations that eliminate or induce 
oscillations. 

More generally, we have enumerated and described the two feedback mechanisms 
that are general to hydrogel coupled chemical reactions: a multiplicative term and 
an inertial feedback. We have shown that these mechanisms can be harnessed to 
design systems that include new positive and negative feedback loops that introduce 
excitable and oscillatory dynamics to an otherwise simple non-oscillatory chemistry. 
We have shown that simple model systems gain control over oscillatory behavior by 
mechanochemically coupling them to the gel. In summary, hydrogel-coupled non-linear 
reactions can be employed to design and engineer complex networks utilizing a full 
toolbox of positive and negative feedbacks. 
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Mechanical pressures 

Chemically cross-linked hydrogels are subject to internal forces arising from both elastic 
stresses and osmotic pressure. The latter is due to both enthalpic solvent-polymer 
interactions as well as the entropically unfavorable displacement of solvent molecules by 
polymer. 

We use the Flory-Huggins model of free energy in a polymer gel to derive the 
internal stress, a, as a function of the polymer volume fraction, 0. For an isotropically 
swollen gel, a is given by the sum of the osmotic pressure and the elastic stress 125; 

^ = -(^ + log(l+x;/;") + Wi„ (7) 

where KT is the thermal energy, cq is the concentration of cross-linking molecules, 0o is 
a reference polymer density, and Vq is volume of one statistical subunit of gel at 00. The 
reference density is that at which the individual sub-units adopt Gaussian conformations 
with minimal elastic energy, y is the polymer-solvent enthalpic interaction term. For a 
PNIPAAm gel Hirotsu [22] gives the enthalpic term as the sum of constant, temperature 
dependent, and density dependent parts: 

X = Xo(r)-Fxi0 (8) 

where xo = 3-42 — 903/r and xi = 0.518. 

The overall enthalpic interaction is augmented by the additional interactions 
between solvent and reactive sidegroups anchored to the gel. In the BZ-gels, these 
are solvent-catalyst and solvent-bipyridine interactions. We assume a general form of 
the interaction that varies with the oxidized fraction of catalyst, fox- 

X Xo(^) T Xl0 T XbpyWat T XrWed T Xox^ox; (9) 

where Ucat is the total number of catalyst molecules per statistical gel subunit, and 
^ox = /ox^cat and Ured = (1 “ /ox)ncat are the number of oxidized and reduced catalyst 
molecules, respectively. 

Under this assumed interaction, we can rewrite the interaction term as = 
(Xo(-C) T ncat(Xred T /ox^Xox) T Xi0) f ; where Xred Xbpy T Xri and Ay ox Xox Xr- 

We next express the oxidized fraction as fox = vVqVo/( ncatf), where v = V/Vq is 
the dimensionless concentration. By substituting for fox, we rewrite the total stress as 

^ = -0 - log(l - 0) - (xo + - X*v(j) + CqUo j , (10) 

where Xq = yo(T) + ^catXred and y* = uoVoAyox- We recover the effective form used in 
previous works caEsi, and expose Equation (9) as their underlying implicit assumption. 
Performing fits to published experimental data !I3. we find approximate parameter 
values of Xred = 12.5 and Ayox = —13.9 
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Mechanical response 


The velocity at the edge of the gel is Wp = ^ where r is the radius of the deformed gel. 
This velocity is proportional to the gradient of the osmotic pressure [IS] : 


3/2 



where Aq is a dimensionless constant describing the mobility of the gel ra Assuming 
a homogeneous gel, we approximate the gradient of the pressure to be Vex ~ cr/r: 



( 12 ) 


For a spherical gel we compute d(/)/dt by applying the relationship (r/R)^ = (^o/^), 
where R is the radius of the gel at (pQ. Manipulation yields 


^ - 4,)a, 


(13) 


where M = 3Xo4>T/R^. 


BZ and Oregonator 

We model the BZ reaction utilizing a modified [IHj and simplified [T] Oregonator. The 
Oregonator captures two key feedback features of the BZ reaction. One feature is the 
autocatalytic production of the species u via the transition of reduced catalyst into the 
oxidized catalyst, v. The second key feature is the inhibition of the reaction via the 
consumption of u by a third reactant, y, which we assume to be in rapid equilibrium. 

The chemical kinetics are described by the two rates dtu = F{u,v,(f)) and dtv = 
eG{u, V, (p), where u and v are reduced forms of U and V, the concentrations of HBr02 
and polymer-bound oxidized catalyst, respectively. The reduced concentrations are 
defined via u = U/Uq and v = V/Vq, where Uq, Vq) fhe characteristic time, and the BZ 
parameters e, /, and q are as defined in [TH] . 

The kinetic functions F and G are 



(14) 


where concentrations are written with respect to total volume, an approximation that 
is valid when the microscopic mixing lengthscale of polymer molecules is much smaller 
than the distance between reacting molecules. 


Oseillatory variant to theoretieal network 

We introduce a variant of the theoretical network presented in the main text. The two 
component model in the main text cannot sustain chemical oscillations but exhibits 
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Figure 7. Integrated trajectory of theoretical model system. The top panel gives 
the concentrations of u (red), v (black), and 2 ; (blue) as a function of unitless time. 
The lower panel gives the polymer volume fraction, (p. The parameters are a = 0.1, 
/3 = 0.2, and M = 1. 


damped oscillations following an initial excitable jump in response to a perturbation 
from the stationary state. This three component variant, however, exhibits stable 
oscillations. Here, both u and v are dissolved in solution and we introduce a new 
polymer-bound reactant, 2 ;. The concentration of bound reactant .2 simply tracks the 
concentration of v: 

F{u, 0) = 1 — u{l — 0)^ 

G{u,v) = auv — I3v (15) 

H{v, z) = /3{v — z) 

where H(v,z) = z. The pressure is a function of the new variable, cr(z). Since .2 is 
enslaved to the auto-catalytic reactant v, increases in v swell the gel. In response, the 
concentrations of the dissolved reactants u and v are increased. In the two variable 
excitable model, the swelling that results from increases in v simultaneously dilutes u, 
creating a drag on the system that damps the oscillations. In this modified model, both 
u and V are concentrated upon swelling whereas .2 is diluted. This system undergoes 
stable oscillations as shown in Figured 
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